home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / ode-initval / rk2imp.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  3.7 KB  |  168 lines
  1. /* ode-initval/rk2imp.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Runge-Kutta 2, Gaussian implicit */
  21.  
  22. /* Author:  G. Jungman
  23.  */
  24. #include <config.h>
  25. #include <stdlib.h>
  26. #include <string.h>
  27. #include <gsl/gsl_math.h>
  28. #include <gsl/gsl_errno.h>
  29. #include <gsl/gsl_odeiv.h>
  30.  
  31. #include "odeiv_util.h"
  32.  
  33. typedef struct
  34. {
  35.   double *knu;
  36.   double *ytmp;
  37. }
  38. rk2imp_state_t;
  39.  
  40. static void *
  41. rk2imp_alloc (size_t dim)
  42. {
  43.   rk2imp_state_t *state = (rk2imp_state_t *) malloc (sizeof (rk2imp_state_t));
  44.  
  45.   if (state == 0)
  46.     {
  47.       GSL_ERROR_NULL ("failed to allocate space for rk2imp_state",
  48.               GSL_ENOMEM);
  49.     }
  50.  
  51.   state->knu = (double *) malloc (dim * sizeof (double));
  52.  
  53.   if (state->knu == 0)
  54.     {
  55.       free (state);
  56.       GSL_ERROR_NULL ("failed to allocate space for knu", GSL_ENOMEM);
  57.     }
  58.  
  59.   state->ytmp = (double *) malloc (dim * sizeof (double));
  60.  
  61.   if (state->ytmp == 0)
  62.     {
  63.       free (state->knu);
  64.       free (state);
  65.       GSL_ERROR_NULL ("failed to allocate space for ytmp", GSL_ENOMEM);
  66.     }
  67.  
  68.   return state;
  69. }
  70.  
  71.  
  72. static int
  73. rk2imp_apply (void *vstate,
  74.           size_t dim,
  75.           double t,
  76.           double h,
  77.           double y[],
  78.           double yerr[],
  79.           const double dydt_in[],
  80.           double dydt_out[], const gsl_odeiv_system * sys)
  81. {
  82.   rk2imp_state_t *state = (rk2imp_state_t *) vstate;
  83.  
  84.   const int iter_steps = 3;
  85.   int status = 0;
  86.   int nu;
  87.   size_t i;
  88.  
  89.   double *const knu = state->knu;
  90.   double *const ytmp = state->ytmp;
  91.  
  92.   /* initialization step */
  93.   if (dydt_in != NULL)
  94.     {
  95.       DBL_MEMCPY (knu, dydt_in, dim);
  96.     }
  97.   else
  98.     {
  99.       int s = GSL_ODEIV_FN_EVAL (sys, t, y, knu);
  100.       GSL_STATUS_UPDATE (&status, s);
  101.     }
  102.  
  103.   /* iterative solution */
  104.   for (nu = 0; nu < iter_steps; nu++)
  105.     {
  106.       for (i = 0; i < dim; i++)
  107.     {
  108.       ytmp[i] = y[i] + 0.5 * h * knu[i];
  109.     }
  110.       {
  111.     int s = GSL_ODEIV_FN_EVAL (sys, t + 0.5 * h, ytmp, knu);
  112.     GSL_STATUS_UPDATE (&status, s);
  113.       }
  114.     }
  115.  
  116.   /* assignment */
  117.   for (i = 0; i < dim; i++)
  118.     {
  119.       y[i] += h * knu[i];
  120.       yerr[i] = h * h * knu[i];
  121.       if (dydt_out != NULL)
  122.     dydt_out[i] = knu[i];
  123.     }
  124.  
  125.   return status;
  126. }
  127.  
  128.  
  129. static int
  130. rk2imp_reset (void *vstate, size_t dim)
  131. {
  132.   rk2imp_state_t *state = (rk2imp_state_t *) vstate;
  133.  
  134.   DBL_ZERO_MEMSET (state->knu, dim);
  135.   DBL_ZERO_MEMSET (state->ytmp, dim);
  136.  
  137.   return GSL_SUCCESS;
  138. }
  139.  
  140. static unsigned int
  141. rk2imp_order (void *vstate)
  142. {
  143.   rk2imp_state_t *state = (rk2imp_state_t *) vstate;
  144.   state = 0; /* prevent warnings about unused parameters */
  145.   return 2;
  146. }
  147.  
  148. static void
  149. rk2imp_free (void *vstate)
  150. {
  151.   rk2imp_state_t *state = (rk2imp_state_t *) vstate;
  152.   free (state->knu);
  153.   free (state->ytmp);
  154.   free (state);
  155. }
  156.  
  157. static const gsl_odeiv_step_type rk2imp_type = { "rk2imp",    /* name */
  158.   1,                /* can use dydt_in */
  159.   0,                /* gives exact dydt_out */
  160.   &rk2imp_alloc,
  161.   &rk2imp_apply,
  162.   &rk2imp_reset,
  163.   &rk2imp_order,
  164.   &rk2imp_free
  165. };
  166.  
  167. const gsl_odeiv_step_type *gsl_odeiv_step_rk2imp = &rk2imp_type;
  168.